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It has been shown recently that predictions from Mode-Coupling Theory for the glass transition of 
hard-spheres become increasingly bad when dimensionality increases, whereas replica theory predicts 
a correct scaling. Nevertheless if one focuses on the regime around the dynamical transition in three 
dimensions, Mode- Coupling results are far more convincing than replica theory predictions. It 
seems thus necessary to reconcile the two theoretic approaches in order to obtain a theory that 
interpolates between low-dimensional, Mode-Coupling results, and "mean-field" results from replica 
theory. Even though quantitative results for the dynamical transition issued from replica theory 
are not accurate in low dimensions, two different approximation schemes -small cage expansion and 
replicated Hyper-Netted-Chain (RHNC)- provide the correct qualitative picture for the transition, 
namely a discontinuous jump of a static order parameter from zero to a finite value. The purpose of 
this work is to develop a systematic expansion around the RHNC result in powers of the static order 
parameter, and to calculate the first correction in this expansion. Interestingly, this correction 
involves the static three-body correlations of the liquid, but more importantly, we show that it 
provides some improvement over the quantitative results of RHNC approximation. 
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I. INTRODUCTION 

A lot of theoretical activity has been devoted to 
the general problem of the glass transition in the last 
decades [1]. Amongst all theories, only two of them 
emerge from a microscopic description, and without as- 
suming the existence of the glass transition itself. The 
first theory is an adaptation of the mode-coupling theory 
of critical dynamics [2] to the glass transition, and was 
devised by Gotze and collaborators [3]. It is nowadays 
simply denoted by Mode-Coupling Theory (MCT). The 
other one is the Random-First-Ordcr Transition theory 
(RFOT) [4], an adaptation to molecular glasses of the 
mean-field replica theory of spin glasses [5, 6]. While 
the development of a first principle theory of dynamics 
is plagued with difficulties (in particular for MCT this is 
mainly due to its lack of self-consistency) [7-11], replica 
theory is a static theory, hence calculations are easier, 
and it has been adapted to the study of the dynamical 
arrest of structural glasses, via three different approxima- 
tion schemes: the small cage expansion [12, 13], the effec- 
tive potential approximation [14, 15] and the replicated 
Hyper- Netted-Chain (RHNC) approximation [16, 17]. 

The first two are formulated in terms of the local cage 
size around a particle. This parameter is expected to 
be infinite in the liquid phase (the particles are free to 
visit all space) but finite and small in the glass phase, 
thus can serve as an order parameter for the glass transi- 
tion. However the large jump of this scalar order param- 
eter (from infinity to small!) seems to prevent one from 
obtaining good quantitative results around the dynami- 
cal glass transition, whereas predictions are more robust 
deep inside the glass phase [12, 13]. The situation is 
reversed for the RHNC approximation: the theory is for- 
mulated in terms of a two-point order parameter (a pair 
correlation function between two replicas, equivalent to 
the non-ergodicity parameter) , that displays a jump from 
zero to a non-zero (but non small) value at the dynamical 
transition. The RHNC approximation correctly predicts, 
in a qualitative way, the phenomenology of the dynam- 
ical glass transition, i.e. the freezing of density fluctua- 
tions due to the appearance of an exponential number of 
metastable states in the system [16, 17]. However, the ne- 
glect of some important corrections (see [15] for a detailed 
discussion) both prevents to obtain accurate quantitative 
results near the glass transition, and gives inconsistent 
results deep inside the glass phase [16]. The main prob- 
lem is that the structural properties of the glass (and in 
particular the non-ergodicity parameter) are found to be 
quantitatively incorrect. For example, the non-ergodicity 
parameter for hard spheres shown in [15] and in Fig. 1 
of this paper poorly compares to the experimental ones 
found for example in [18]. 

The difficulties of the replica method in treating the 
dynamical glass transition, compared to the quantitative 
success of MCT in the same regime [19], must not hide 
its successes: it is able, within a purely static framework, 
and starting from the Hamiltonian of the system, to pre- 



dict without assuming it the appearance of a large num- 
ber of metastable states, to quantify this number through 
the complexity, and to make predictions for the critical 
properties of the long-time dynamics [40-42]. It is thus 
of primary importance to continue developing the RHNC 
approximation scheme in order to probe the limits of this 
replica (static) approach, and in order to set the stage for 
an hypothetical theory that would be able to reconcile dy- 
namic, MCT-likc computations, and static, replica com- 
putations. Interestingly, a first bridge was made between 
the two [20] by looking at the dynamics in a replicated 
liquid, and treating these dynamics in a MCT fashion. 
Unfortunately, this approach suffers from the same lack 
of internal self-consistency as the original MCT compu- 
tation, and it only succeeds in reinforcing the general be- 
lief that MCT and replica theory can indeed be unified 
in a coherent scheme (still neglecting activated events), 
as was done in the context of spin glasses [21]. 

In this paper we perform a first step in the general 
task of making RHNC an efficient quantitative theory of 
the dynamical transition, by setting up a perturbativc 
expansion in powers of the static order parameter used 
in replica theory to characterize the glass phase. This or- 
der parameter is proportional to the non-ergodicity fac- 
tor, and our expansion thus corresponds to the "weak 
glass" expansion sought for in [22]. We calculate the first 
correction to RHNC in this expansion, and show that 
it already leads to some improvement of the quantita- 
tive results, while the qualitative picture remains stable. 
Interestingly, we show that this new term, and the subse- 
quent ones, involve the three-body and higher order cor- 
relations of the liquid, which currently receive renewed 
interest: even though two-body static correlations of the 
liquid are blind to the presence of the glass transition, 
it is possible that three-body functions are more sensi- 
tive to it [23] . Our work paves the way to more powerful 
re-summations, maybe allowing to unify the three ver- 
sions of replica theory for structural glasses (small cage, 
effective potential, and RHNC) in a unique scheme, and 
possibly obtain a quantitatively competitive theory of the 
dynamical transition. 

The paper is organized as follows. In Section II we 
briefly review the replica interpretation of the glass tran- 
sition, focusing on the two main order parameters, static 
and dynamic, currently used to describe it in a theoreti- 
cal way. In Section III we introduce the replicated liquid 
theory that will allow us to perform approximations on 
the static order parameter, and briefly review its simplest 
implementation, the replicated Hypcr-Ncttcd-Chain ap- 
proximation [16]. We then explain how to proceed from 
this approximation to perform an expansion in powers 
of the static order parameter. In Section IV, the first 
correction to the replicated Hyper-Netted-Chain approx- 
imation is obtained, and shown to involve three-body 
correlation functions of the liquid. This section is mostly 
technical and can be safely jumped by the reader uninter- 
ested in details of the computation. A numerical solving 
of this improved Hypcr-Ncttcd-Chain approximation is 
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presented in Section V, and we conclude in Section VII. 



II. STATIC ORDER PARAMETER FOR 
GLASSES 

The glass transition is a priori defined in a dynamical 
way as the temperature below which (or density above 
which) density fluctuations become frozen in the liquid. 
However, replica computations on mean-field spin glasses 
have given an insight into the thermodynamics of glasses, 
and provided us with a static order parameter for the 
glass transition [6]. In this section we recall the proce- 
dure needed to define and calculate this order parameter, 
and show its relation with the dynamical order parame- 
ter traditionally used in glass physics, namely the non- 
ergodicity parameter. 



A. Order parameter with replicas 

In order to detect the glass transition from a static ob- 
servable, one introduces identical copies of the system of 
interest . Consider an integer number m of copies of the 
system, i.e. a liquid of m x N particles. The copies are 
indexed by alphabetical characters a, 6, . . ., and the posi- 
tion of particle i in copy a will be denoted by xf. Within 
each copy, all the particles interact with a given pair po- 
tential, but each particle in a given copy is attracted to 
all other particles of all other copies via an attractive 
pair potential of infinitesimal amplitude e. Consider for 
example copy a and copy b, with a ^ b. In analogy with 
usual liquid theory (see [24] and Eq. (15) in the follow- 
ing), each of them has a microscopic density defined by: 



N 



(1) 



• if p < pd and letting e go to zero, the two copies 
de-correlate, i.e. h a b = 0, and the system is not 
trapped in a metastable state. In this case the limit 
e — > and N —> oo commute. 

• if p > pd, the two copies are trapped into the same 
metastable state by their mutual attraction e. Let- 
ting £ go to zero after N — > oo, the two copies will 
stay correlated and h a b ^ 0. Note that obviously 
the copies will de-correlate if e is set to zero before 
N — >• oo, hence the two limits do not commute in 
the glass phase. 

Of course this transition is observed only if metastable 
states are present in the system, that are able to constrain 
copies to stay in the same state, in the thermodynamic 
limit and in the long-time limit. As we will see in the fol- 
lowing, this assumption can be checked self-consistently 
within the theory. In order to recover the equilibrium 
properties of the original, non-replicated system, one has 
to make an analytic continuation to non-integer values of 
m, the number of replicas, and take the limit m — > 1 at 
the end of the calculations [16, 26]. 

The general procedure is then the following: 

• Compute the equilibrium properties of a m-times 
replicated liquid, with small attraction e between 
different copies 

• Take the thermodynamic limit, and compute the 
desired quantities as functions of m 

• Send the attraction e to zero 

• Take the limit m — > 1 to recover the original system 

We show in the following that calculating the static 
order parameter defined above with this prescription 
amounts to evaluating the non-ergodicity parameter of 
the glass. 



Consider now the generalization of the pair correlation 
function, as defined in [24] and Eq. (16) below, to a repli- 
cated system: 



, / \ (Pa(x)pb(y)) . , , , 
h ab (x,y) = 3 1 for a b , 

P 



(2) 



Note that, for a ^ b, the coincident point term is absent, 
because particle i of copy a interacts with all particles 
of copy b, including particle i, whereas within copy a, 
particle i does not interact with itself. 

We want to study what happens in the thermodynamic 
limit N — > oo when the coupling potential strength e is 
vanishingly small. As usual, a phase transition will be 
signaled by the fact that the limits N —> oo and £ — > 
do not commute [25-27]. More precisely, if we suppose 
that a dynamical transition occurs at p^, with a mean- 
field phenomenology, we will have that: 



B. Link with the dynamic order parameter 

In the mean-field replica interpretation, the partition 
function of the system is supposed to be separated into 
many pure states, that correspond to minima of the free- 
energy. Because of this separation, averages can be sep- 
arated in two operations: first an average inside a state 
a, denoted by (•) Q , and then an average over all states 
denoted by • . Calculating h a b(x,y) e , the correlation 
between copy a and b (with a ^ 6), we take into account 
that the attractive coupling between a and b will force 
them into the same state a leading to: 



hab{x,y) e 



(p{x)p(y)) c 
P 2 



1 



(3) 



Letting the interaction go to zero will allow the two repli- 
cas to dc-corrclatc inside the state, leading to: 



hM^) t = tf(4 f (4 -i 



(4) 
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Note that the average density inside a state need not 
be constant, because translational invariance is restored 
only after summation over all the states. 

Independently, the time-dependent density-density 
correlation that is the focus of MCT, and more gener- 
ally dynamical calculations, is defined by: 



In this paper, we will consider for concrcteness a sys- 
tem composed of N spheres of diameter D interacting 
via a hard-sphere pair potential v: 



v(x,y) = 



oo if \x — y\ < D , 
otherwise . 



(12) 



F(x,y,t) = ((p(x,t) - p)(p(y,0) - p)) , 



(5) 



where at t = the system is taken at equilibrium. Note 
that the initial time value of this function is related to 
the inverse Fourier transform of the structure factor: 



F{x, y, 0) = p 2 h(x, y)+5(x~ y)p . 



(6) 



Above the glass transition, and in the long time limit, 
the system eventually gets stuck in a mctastablc state a, 
and we have: 



F{x,y,t) -> (p(x)p(y)) a - p 2 

t—¥O0 



(7) 



But the system is at least able to de-correlate inside the 
state, giving: 



F(x,y,t) -> (p(x)} a (p(y)) a -p z = p 2 h ab (x,y) . (8) 

Now taking the Fourier transform of equations (6) and (8) 
we get: 



F(k,t) -> P 2 h ab {k) , 

t— >oo 

F(k,0)=pS{k) , 



(9) 



where S(k) is the structure factor defined in [24] and 
Eq. (17) below. Finally the non-ergodicity factor is tra- 
ditionally defined as: 



f(k) = lim 



F(k,t) 



*->oo F(k,0) 



so that we get [15, 20]: 



f fc ) = c/m Wlth a r b ■ 
S(k) 



(10) 



(11) 



We thus find that, in the replica interpretation, the dy- 
namical and static order parameter are the same physical 
observable. This parallel is exploited in order to compute 
long-time properties of the dynamics from replica calcu- 
lations, for example in [40] or [41, 42]. 



III. EXPANSION IN POWERS 
OF THE ORDER PARAMETER 

Before trying to obtain a good theory for our replicated 
system and studying the glass phase, we must introduce 
the liquid theory tools used to describe the liquid phase, 
in order to fix notations and for the sake of completeness. 



Note that we use here the hard-sphere potential mainly 
for the sake of simplicity, because it has only one param- 
eter, the density, and because its liquid properties are 
well studied, both numerically and theoretically. How- 
ever the general computation scheme we develop here is 
not restricted to this potential and can be applied to any 
pairwisc additive potential. 

A microscopic configuration a of the system spec- 
ifies all the positions Xi of the particles, with 
i G [1;-^]- The Hamiltonian of the system is 
then H(a) = X^<j v i x i — x j)- F° r anv inverse tempera- 
ture j3 = 1/ksT, the Boltzmann weight of such a config- 
uration is e~P H ( a \ which can be either if there is any 
overlap between spheres, or 1 if no spheres overlap. The 
temperature thus plays no role here. To stay in contact 
with standard liquid theory notations, we will nonethe- 
less keep it explicit in this subsection. We fix units so 
that D = 1 and &b = 1. 

For a given chemical potential p and an eventual exter- 
nal potential ^, the grand canonical partition function of 
the system is defined as: 



Z liq = Tr e-^«^})+£^) , 



(13) 



where the grand-canonical trace is defined 

as Tr • = J2n=o W. I Hi=i dXi *' an< ^ tne general- 
ized chemical potential is defined as v(x) = (3p — f3^{x). 
The grand-canonical average of an observable O is 
defined by: 

(0({xi})) nq = ^Tr 0({*,})e-/^>H^K*O (u) 

We define the microscopic density p and the average den- 
sity p as: 



N 



p{x) = (p(x)) n . 



(15) 



We limit ourselves to translationally invariant systems, 
so that the average density is constant in the system and 
p(x) — p in all points. The pair correlation h is linked 
to the normalized second cumulant of the microscopic 
density, apart from a coincident point term: 



h(x,y) = -j (y^^TSix - Xi)8(y - Xj) \ -1 



=1 fri 



liq 



(p{x)p(y))n q i 1., A 
5 ^ - 1 - -S(x - y) . 



(16) 
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Subtracting 1 gives a function that decays to zero at 
large interparticle distance, and thus allows to compute 
its Fourier transform. In a translational invariant sys- 
tem, h(x,y) depends only on \x — y\, and we define the 
structure factor S. related to its Fourier transform as: 



S(k) = 1 + ph{k) 



(17) 



h is only the second function of a whole hierarchy of 
cumulants of the microscopic density. In the following 
we will need the next cumulant in this hierarchy: 

W£\x, y, z) = (p(x)p(y)p(z)) + 2p(x)p(y)p(z) 



p(x) (p(y)p(z)) 



{2 permutations} . 

(18) 



Usual liquid theory approximations start form the Legen- 
dre transform of In Z with respect to the chemical poten- 
tial, using the fact that the chemical potential is coupled 
to the microscopic density. Thus we define: 



r iiqW : 

is* [p] such that 



p(x)is*[p](x) -\nZ iici [v*[p}] , 
<SlnZ liq [i/] 

= P[x) ■ 



5is a (x) 



(19) 



v *[p) 



From the virial expansion [43], one can show that this 
functional is composed of an ideal gas term plus an 
excess term: 



r iiq[p] =r id [p] + r ex [p] , 



p{x) [lnp(x) - 1] 



(20) 



From this functional, one can define a hierarchy of corre- 
lation functions (equivalent to the vertex functionals in 
field theory) called the direct correlation functions [24] : 



( n ) i \ 

C liq i X U--- ,XN) 



s n (r liq M-r id [ P ]) 

Sp{xi) ■ --SpixN) 



(21) 



By the properties of the Legcndrc transform (the calcula- 
tion is shown formally in Eq. (57) and can be straightfor- 
wardly applied to the case of liquid theory) , the second 

(2) 

order direct correlation function c liq , that will be denoted 
by c in the following, is the inverse function of h, which 
is expressed by the Ornstein-Zcrnikc equation: 



(l + ph(k))(l-pc(k)) = l 



(22) 



( 2) 

where c(k) is the Fourier transform of c liq (r) as defined 
in Eq. (21). We can rewrite this equation in a more usual 
field-theoretic notation: 



G Uq (x, z)T ( ^ q {z, y) = 8(x, y) 



(23) 



where wc have defined the propagator Gu q (x,y) ~ 
p(x)8(x 1 y) + p(x) p{y)h{x , y) and the two-point vertex 



functional of liquid theory T^(x,y) = -^8(x,y) - 
c(x,y). Differentiating this equation with respect to the 
density we obtain the third order Ornstcin-Zcrnicke equa- 
tion: 

,(3) 



{x, x')T ( ul (v, i/')rg (*, w$ (*', y ', z') 



(24) 



(3) 

where T,- is related to the third order direct correlation 



liq 



function via Eq. (21): 



Hiq( x >y> z ) = -^^ s (x,y)S{x, z) ~ c\'^{x,y, z) . (25) 

We will make repeated use of relations (23-25) in the 
next section. 
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(3), 



A. Replicated liquid theory 

We can now formulate in the same way the partition 
function of our replicated system, defining: 



Z m = Tr e 



(26) 



where v a b is equal to v for a — b, and v a b is a small attrac- 
tive coupling when a b. The prime on the summation 
sign ^2 means that when a = 6, the summation must 
exclude the case i — j , and the trace operation is now de- 
fined as: Tr • = ^^=0 M™" / UT=i UZ*i dx t •• For thc 
sake of simplicity, u in the following will denote the set 
of m chemical potentials {va}a=i...m- an d v will denote 



the set of m pair potentials {v a b} a . 



6=1. 



Equivalently 



the family of m average densities {p a }a=i...m defined in 
Eq. (1) will be denoted by p and the family of m 2 corre- 
lation functions {h a b} a ,b=i...m will be denoted by h. 

As for the non replicated liquid case, one can define the 
Legendre transform of the free-energy with respect to v, 
to obtain an m-dependent Legendre transform Ti^ m : 



r i,m[p] = / P^ X ) V 1[P\{ X ) - faZ m [v*] 



SlnZ m [u] 



(27) 



5v a {x) 



Pa{x) . 



However, our goal here is to focus on the static order 
parameter h a b defined earlier. It is thus preferable to 
obtain a theory that explicitly depends on it in order to 
control approximations in terms of this quantity. The key 
point is to notice that the pair potentials v a b are coupled 

(2) 

to the two-point densities of the liquid p ab , defined by: 



Pab'( X >y) = Pa{x)pb{y) - Pa(x)S a b5(x i y) , 

p ( ab( x >y) = (Pab&y) 



(28) 
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This quantity is trivially related to h a b by p^ b \x,y) = 
Pa{x)Pb{y) [1 + h a b(x,y)], and we define it for conve- 
nience. Indeed we remark that: 



5 In Z r , 



8(-0Vab(x,y)) 
S/3v ab (x,y) 5 



1 ( 2 )r ^ 



(29) 



This relation explains why it will be easier to manipulate 
the properties of the Legendre transform that we define 
below in terms of p^ than in terms of h, avoiding unnec- 
essary density factors in the calculations. Also for later 
convenience, we define the propagator and two-point ver- 
tex function: 

Gab{x, y) = p a {x)5 a bS(x, y) + Pa{x)pb{y)h a b{x, y) , 



1 



Pa(x) 



5 a bS(x,y) - c a b{x,y) 



(30) 



where c ab is the generalization to mixtures of the direct 
correlation function. Again we have the OZ relation, i.e. 
that Tt T^c (x, z)G c b{z,y) = 8 a b5(x,y), or equivalently 
h a b(x, y) = c ab (x, y) - ^ c J z h ac (x, z)p c (z)c cb (z, y). 

Defining T m the Legendre transform of Fi. m with re- 
spect to /Jv we get: 



r m [p,h] = --^ / p^ b \x,y)f3v* b (x,y) 

L a.b Jx >y 
+ r 1 , m [ P ,v*\ , 



(31) 



<n\ m [p,v] 



5f3v ab (x,y) 



^Pab( x >y) > 



where we have written T m as a function of h since p^ 
and h are simply related. 

The functional T m , when evaluated at the true correla- 
tion function, coincides with Li iin , but now a stationary 
principle allows to derive self-consistent equations for the 



correlation function, because as a consequence of Legen- 
dre transforms properties we have: 



ST m [p,h] 

Sh ab (x,y) 



Pa(x)p b (y) 



sr m [p, h] 



(32) 



■^Pa{x)p b {y)l3v* ab [h\{x,y) 



We know that the values v* b of the pair potentials that 
lead to the real correlation function are v ab , so that the 
real correlation functions h ab are solutions to the self- 
consistent equation: 



£r m [ P ,h] 

Sh ab (x,y) 



1 



Pa(x)p b {y)f3v ab (x,y) 



(33) 



Now performing an approximation for L m and using 
Eq. (33) leads to an approximate value of h a b, and eval- 
uating the approximated L m at this value of h ab lead to 
an approximate free-energy that is consistent with the 
approximate h ab obtained. Finally, as described in Sec- 
tion II, we set the inter-replica potential to zero, which 
distinguishes between the case a = b, where v ab = v, 
and the case a ^ b, where v ab is now set to zero. Note 
that we set the inter-replica potential to zero only after 
taking the Legendre transform. This ensures that the 
thermodynamic limit has been taken, and will give rise 
to non trivial inter-replica correlations. We look for per- 
sisting inter-replica correlations, which will indicate the 
glass phase [16]. 



B. Morita & Hiroike functional 

Morita & Hiroike [44] showed that the functional F m 
can be written as: 

T m [p, h] = T Id [p, h] + T mng [p, h] + F 2P/ [p, h] , (34) 



J 



where: 



rjd[p,h] = y~] f p a (x) [lnp a (x) - 1] + ^ y~] I p a (x)pb{y) [(1 + h ab (x,y))ln(l + h ab (x,y)) - h ab (x,y)} , 

a Jx a,b Jx 'V 



1 (—1)™ 

T Ring [p,h] = 2 X! ~ Tl " Pai(xi)h ai a 2 (xi,X 2 ) ■ ■ ■ Pa n (x n )ha nai (x ni Xl) , 



(35) 



n>3 



and L2p/[p, h] is the sum of all 2PI diagrams, which are 
usual Mayer diagrams composed of black nodes p a {x) 
and links h ab (x, y), such that when two links are removed 
from the diagrams, it does not disconnect in two separate 
parts, r Ring is the sum of all ring diagrams. The Hyper- 
Netted-Chain (HNC) approximation [45], and other per- 



formant approximations of liquid theory start from this 
exact functional. This result is exactly what we wanted: 
we have now an explicit functional of our static order 
parameter. Before turning to the expansion of this quan- 
tity in powers of the order parameter, we first review the 
results obtained within RHNC approximation, starting 
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from the above Gibbs free-energy, since it is at the core m — > 1. In order to do that, we will assume replica 
of our expansion. symmetry, i.e. that: 



C. RHNC approximation 

RHNC amounts to discard all the 2PI diagrams, giv- 
ing an analytic expression for the free-energy as a func- 
tional of h which is just the sum of the two terms in 
Eq. (35). With this approximate functional, that we will 
denote by T^ NC , we obtain self consistent equations for 
the h ab by making use of Eq. (33). The sum of ring 
diagrams, when differentiated with respect to h ab , i.e. 
when cutting a link, gives a sum of open chains, which 
is equal to h — c. Indeed, rewriting the generalization of 
the Ornstcin-Zcrnikc equation (22) to multicomponcnt 
mixtures in direct space we have: 

Cab(x,y) = h ab (x,y) - h ac (x, z)p c (z)c cb (z,y) , 

c Jz 

(36) 

which can be solved iteratively with respect to h, to give: 
c ab (x,y) = h ab (x,y) 

OO 

- }X-l) n ha ai (x, X 1 )p ai (Xl) ■ ■ ■ Pa n (Xn)ha n b(x n ,b) , 
n=l 

(37) 

where summation over repeated indices and integration 
over repeated positions has been assumed. Comparison 
with Eq. (35) shows that h ab — c ab is indeed exactly the 
derivative of Tm ng with respect to h ab : 



8T Rlng [p, h] 
8h ab (x 7 y) 



1 



Pa{x)pb(y) [h a b{x, y) - c ab (x, y)} 



(38) 

Now using the variational principle stated in Eq. (33) 
along with the prescription described in Section II, we 
get a set of m 2 self consistent equations: 

ln(l + h aa (x, y)) = -/3v(x, y) + h aa (x, y) - c aa {x, y) , 
ln(l + h ab (x, y)) = h ab (x, y) - c ab (x, y) for a ^ b . 

(39) 

These equations must be supplemented with the repli- 
cated Orstein-Zernike equations (36) that read in Fourier 
space: 

Caa(k) =h aa (k) - ph aa (k)c aa (k) - p^h ab (k)c ab (k) , 

b^a 

c a b(k) =h ab (k) - ph aa (k)c ab (k) - ph ab {k)c bb (k) 
~P ^ h ac (k)c cb (k) . 



h aa {x,y) = h{x,y) Va , 
h ab (x, y) = h(x, y) Va ^ b . 



(41) 



This allows to perform the summations over replica in- 
dices in Eqs. (39-40) and perform the m — > 1 limit, to 
obtain two sets of equations. In the first set, we see that 
the functions h and c decouple from h and c to give the 
usual liquid theory HNC approximation and Ornstcin- 
Zernike equation: 

ln(l + h(x, y)) = -/3v(x, y) + h(x, y) - c(x, y) , 
c(k) = h(k) - ph(k)c(k) . 

_(42) 

The second set of equations specify the functions h and 
c, given h and c: 



ln(l + h(x, y)) = h(x, y) - c{x, y) , 

c(fc) = h(k) - ph(k)c(k) - ph{k)c{k) + ph(k)c(k) 



We can rewrite the second of these equations as: 



c(k) 



(1 - pc(k))h(k) 
1 + ph(k) - ph(k) 



(43) 



(44) 



and use the definition of the structure factor Eq. (17), the 
Ornstcin-Zcrnike equation Eq. (22) and the link between 
h a b for a ^ b with the non-crgodicity factor Eq. (11) to 
obtain: 



i 



/(*) 



P S{k) 1 - f(k) 



(45) 



Since c and h are both related to /, we see that Eq. (43) 
is a self-consistent equation on the non-ergodicity param- 
eter, which reads 



1 



/(*) 



P S{k) 1 - f(k) 



F[h - \n{l + h)](k) , (46) 



where J- denotes a Fourier transform and h is expressed 
in terms of / by Eq. (11). If we were to make an ex- 
pansion of the r.h.s. of this equation in powers of f(k) 
(hence of h), we would get, at lowest order: 



f(k) S(k) 



(40) 



1 - f(k) 2p 



S(q)S(k-q)f(q)f(k-q) . (47) 



In order to follow the prescription of Section II, we must 
make the free energy analytic in to, then take the limit 



This form of self consistent equation is very reminiscent 
of the Mode-Coupling result that read [28]: 
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f(k) S(k) 



Muc(k,q)S(q)S(k-q)f(q)f(k-q) , 



1 - f(k) 2p J' 

2 2 

Muc(k, q) = (ft ■ «c(«) + k • (ft - g)c(ft - q) + pk 2 c^(k, -q)J 

= [1 " ^(9)1 + fc " ( ^~ g) [1 - Mfc - ?)] - [l + p 2 cg(fc, - ? ) 



A- 2 

ft 2 



ft 2 

^Vffife) + ^t^prg(ft - ff) + Ag(*, -*) 



(48) 



ft 2 



Note here the presence of the three-body direct corre- 
lation function of the liquid, even though it is usually 
neglected, since it has been shown to be negligible with 
respect to the other, two-body term [29], except for spe- 
cial cases [30] . However the RHNC result if expanded at 
this order in h a b is trivial (in the sense that the kernel 
A4 = 1), and reminiscent of the result obtained by [9] 
by a dynamical field-theory calculation. Whether this is 
a coincidence or not is an open question at this stage. 
Nevertheless, we performed this expansion of the loga- 
rithmic term in Eq. (43) only to show the striking simi- 
larities that exist between replica calculations and usual 
Mode-Coupling ones. 

If we do not perform this truncation, and solve the 
RHNC approximation Eq. (43) combined with the HNC 
approximation for the liquid part Eq. (42), we already 
find, without assuming it, the existence of a dynamical 
transition [16]: for hard spheres, for densities lower than 
tfd ~ 0.599 (where the packing fraction tp = pit/ 6), the 
solution to Eq. (43) is always h a b = 0, i.e. a liquid phase, 
whereas h a b discontinuously jumps to a non zero value 
for tp > tfd, indicating the glass transition. This results 
is quite good, because the dynamical transition for three- 
dimensional hard spheres is estimated to be around ipd « 
0.57. Note that Mode-Coupling Theory instead strongly 
underestimates the transition [19]. We show in Fig. 1 the 
resulting non-ergodicity factor [15] for a packing fraction 
(p = 0.6, obtained by solving the set of equations (43). 
This is the crucial problem of RHNC: f(k) is quite far 
from the numerical results (see [18, 19]) which are instead 
well captured by the Mode-Coupling theory. 

The goal of this paper is to demonstrate that the order 
parameter h a b with a ^ b can be used as an organizing 
device for the theory in order to gradually incorporate 
higher-order correlations of the liquid into the replica re- 
sult. Of course, we see in Fig. 1 that the order parameter 
is not a small quantity, and thus an expansion in pow- 
ers of h a b is not a priori justified. Note however that 
RHNC already re-sums an infinite number of diagrams 
containing arbitrary numbers of h a b links, which maybe 
explains its ability to predict a transition towards a non- 
small value of the order parameter. Our purpose here is 
to build from RHNC and incorporate more diagrams, and 
we will see that even keeping the lowest order correction 
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FIG. 1: Filled squares: Non-ergodicity factor of hard spheres 
at packing fraction 0.599670 (i.e. at the transition), ob- 
tain from the replicated Hyper-Netted-Chain approximation 
Eq. (43), combined with the liquid theory HNC approxima- 
tion Eq. (42). Open circles: the corresponding order param- 
eter h(k) 

already provides some improvement over RHNC. 

D. Improvements over the liquid quantities 

Before turning to the study of the expansion in powers 
of the order parameter, a remark is in order. We have 
seen that in the m = 1 limit, the liquid quantities de- 
couple from the inter-replica correlations: Eq. (42) is in 
fact independent from the off-diagonal correlations. This 
allows us to use any liquid theory approximation to eval- 
uate liquid quantities that appear in our equations. For- 
mally, this is justified by writing the Gibbs free-energy T 
in Eq. (34) as: 

r m [ / ,,h] = rf c [ / ,,h] + r^[p,h]+r^i s [p,h] , (49) 

where T 2 ^ 1 is the sum of all 2PI diagrams that do not 
contain any h a b, a ^ b links. We can now use the varia- 
tional principles, which will give the full liquid correlation 
function for the a = b components. For a given approxi- 
mation of r 2 ^ s , we will obtain a self consistent equation 
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FIG. 2: Non-ergodicity factor as a function of the wave vector 
at packing fraction 0.6, for the standard replicated HNC cal- 
culation, and from the combination of PY approximation for 
the diagonal part and RHNC for the off-diagonal correlation. 

for the a ^ b components. Neglecting altogether T^^, 
we recover the set of equations Eqs.(43), in which c and h 
are now the full liquid correlation functions. Of course, 
inside the glass phase, the liquid quantities cannot be 
obtained numerically or experimentally, so that such a 
full re-summation is useless. Instead, we need to have an 
approximation that can be extrapolated from the liquid 
phase. For example we can choose the PY approximation 
for the liquid quantities. Choosing to work at a packing 
fraction (defined as (p — Trp/6) tp = 0.6, which is ex- 

I 



where at ^ 6; Vi, and we have underlined the fact that 
the zeroth order term is absent since the full functional 
must vanish in the liquid phase, where T^f ass — by 
construction. 

Now consider the "glassy" 2PI diagrams in Eq. (49). 
We can show that these diagrams must contain at least 
three h lines: this idea was already used in the "weak 
glass" expansion of [22, Appendix A3]. Indeed, a h link 
joins two nodes that have different replica indices, say a 
and b. All the nodes connected to the a node by a path 
of h links must also have replica index a, and the same 
applies for the b node. Thus all h links are nodal links: 
they separate the diagram in two parts, each of it has a 
different replica index. If a 2PI diagram would contain 
one or two h links, then differentiating once or twice with 
respect to h would cut the diagram in two, which is in 
contradiction with the fact that the diagram is 2PI. Thus 
we proved that all 2PI diagrams contain either zero or 



pected to be above the dynamical transition, we solved 
the RHNC equation for the off-diagonal part, Eq. (43) 
with the PY direct correlation function as an input for 
the diagonal part, and found that this brings about a 
little improvement over standard RHNC results. It is 
known that the PY approximation gives a less important 
underestimation of the peak of 5(fc), which is the main 
ingredient that leads to the glass transition. In standard 
RHNC as in the PY version of it, we find a transition 
from a liquid state at (p < tpj to a glass state at ip > (pd, 
where the self-consistent equation on h admits a non- 
zero solution. The value of the critical density is shifted 
downwards from ip c w 0.599 to 0.591 by the use of PY 
approximation, which is an improvement, even if mod- 
est. We show the resulting non-ergodicity factor in the 
two sets of approximations in Fig. 2. Even though use of 
the PY approximation gives a slightly larger non ergodic- 
ity factor, the result is still very small when compared to 
simulation and experimental data, where / is much closer 
to 1 at small wave vectors. We conclude that the source 
of the problem in the static approach is not the diagonal 
part of the liquid. We should therefore seek for a way of 
improving the equation for the off-diagonal correlation, 
Eq. (43). 

E. Systematic expansion in powers of h 

To go further than the replicated Hyper-Netted-Chain 
approximation, we perform an expansion of ^ g fJ ss in 

powers of h: 



(50) 

I 

three or more h links. Moreover, a diagram that contains 
three h links can have at most six parts composed of h 
links and p nodes that all have the same replica index. A 
little reflexion shows that the only possibility to construct 
such diagram and make it 2PI is the one pictured in Fig. 3. 



From this analysis, we showed that the non HNC terms 
in the free-energy Eq. (34) can be organized as a series 
in powers of h, the first term of which are 0(h 3 ). This 
last observation is crucial, since it will allow us to make 
the computation tractable. 



In the following we will need to distinguish the deriva- 
tives with respect to the density and with respect to the 
correlation functions, as well as derivatives with respect 
to chemical potentials and with respect to pair potentials. 



°° 1 8 n T 2PI J^kj 

a gla SS [pM = y^Tr ' — — g ' ass ' rh aibl (xi,yi) ■ • ■ h anbn (x n ,y n ) , 
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FIG. 3: Diagrams that contribute to the free-energy at order 
h Li . A wiggly line joining two replica indices a and b is a 
h a t, with a ^ b function, a black dot attached to a zone with 
replica index a is an integration point weighted by a density 
factor p a - 



We define: 



sr m [ P ,h] 

5pa{x\) 

ST m [p,h] 

6h ab (xi,yi) 

<5 2 F m [ P ,h] 

8p a {xi)8h cd {x2,y2) 



(51) 



and so on. We will also need to define derivatives of 
In Z m with respect to the chemical potentials and pair 
potentials: 



SlnZ m [v, v] 



w. 



Hi, l ' 1) (^i,yi) 



5v a {xi) 
8lnZ m [v, v] 



6(-/toab(zi,2/i)) 
5 2 lnZ m Kv] 



<5i/ a (xi)(5(-/fr; C( i(x2, y 2 )) 



(52) 



and so on. Because of our analysis, we see that the n = 1 
and n = 2 terms in Eq. (50) are necessary zero. In the 
following, we will calculate the first non-zero term in this 
expansion, the third-order term. In practice, we found 
it easier to calculate the third derivative of the total free 
energy, and substract from it the third order term of the 
RHNC free-energy, i.e. we calculated, with the above 
notations: 



-p(°.3) i \ 
K b ' iC ' dief {xi,yi-,x 2 ,y2]X3,y3) 



6 3 



^HNC 



(53) 



5h ab (xi,yi)8h c d(x 2 ,y2)Sh e f(x 3 , y 3 ) 



with a ^ b, c =^ d and e ^ /. The third order term 
of the RHNC free-energy will be simply calculated from 
Eqs. (34-35). 



IV. TREATMENT AT THIRD ORDER IN THE 
ORDER PARAMETER 

In order to make progress, we must find a way to eval- 
uate the third derivative of the free-energy with respect 
to the correlation function. We can make use of the prop- 
erties of the Legendre transform to do so. The only diffi- 
culty is that we performed two Legendre transformations, 
with respect to different objects, which renders difficult 
the use of properties such as the one written in Eq. (24) . 



A. Third-order Ornstein-Zernicke relation 

Considering, for simplicity, a discretized version of our 
theory, we have that p is a m by M matrix, where M is 
the number of points of the underlying lattice, and h is 
a m by m by M by M object, and the same applies to 
v and v. We can write the two pairs p, p^ and v, v in 
two (big) vectors: 



* - [p, ^P (2) 

* = (U, -p\) 



(54) 



with the convention that if i > mM an index i of the 
vector * or $ must be understood as a group of two 
spatial coordinates (xi, yi) and two replica indices (a^, bi), 
but if i < mN, it must be understood as one spatial 
coordinate and one replica index. The double Legendre 
transform F m can then be written as: 



r m [*] =Tr $**-lnZ m [**] , 
81nZ m [<f>] 



with ** such that 



<5* 



(55) 



= * 



We have by definition of the Legendre transform: 

sr m m 



(56) 



We can perform a functional derivative of this equation 
with respect to *2 to get: 



S 2 lnZ m {$] 



<5*i(5$2 



5*1(5*2 



^ <5i2 = Tr 



<5 2 r m [*] <5 2 lnZ m [$] 



<5*i(5* 4 (5$ 4 (5$2 



(57) 



Computing a third derivative with respect to *3, we get 
formally: 



= Tr 



Tr 



/ <5 3 F m [*] <5 2 lnZ m [<i>] 



\ (5*1(5*3(5*4 <5*2f>*4 

<5 2 F m [*]<5 2 r m [*] <5 3 lnZ m [$] 



<5*k5*4 (5*4(5*5 (5*2<5<I'4(5<J) E 



(58) 
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And multiplying through by a second derivative of T r 
and using Eq. (57), we get: 

5 3 r m [*] 



5*i5* 2 5* 3 

/5 2 I\„[*] 5 2 r m [*] 5 2 r m [*] 5 3 lnZ m [*] 



Tr 



\ J^id*!' 5* 2 5* 2 ' 5* 3 5* 3 ' 5*i'5* 2 '5*3 



(59) 



This is a standard equation that can be found in any text- 
book on field theory, for example in [31], that expresses 



the relation between vertex functions and propagators of 
the theory. It coincides with Eq. (24) if the vectors * 
and * contain only the first mM components. 

We are interested in the continuum space limit of this 
expression, evaluated in the liquid, and with 1,2 and 3 
that are all greater than mM. With explicit spatial coor- 
dinates and replica indices and implicit summation and 
integration over repeated indices and variables, and after 
replacing derivatives with respect to p^ with derivatives 
with respect to h, we get: 



= (*i,vi\ z'OriS (X2,y2-, 4)r2# (* 3 , y 3 ; x' 3 )w^ e , K; 4; 4) 

16 ^(0,2) / , M-r(M)/ n„,(2,l) ' / / ' I \ 

Kb,a'b'( X ^y^ X 2^ 3 )T y cd J {X2,y2]X 2 Fef,e< 0=3, lfi»5 x sWa'¥ ,c' ,e' { X vVv X 2> x 3> 

- two permutations {(a, b); (xi,yi)} <H> {(c, d); (x 2 , 2/2)} <-> {(e, /); (0:3, ys)} 

32 t-,(1,1)/ /\t^( > 2 ) / ' t\-n(°' 2 ) I 1 /\ii/(1-2) ' / ' / ' /\ 

- Tio V ab,a\ X ^y^ X l) V cdJ>d>\ X 2^ 

- two permutations {(a, b); (#1,2/1)} o {(c, d); (x 2 , y 2 )} O {(e,/); (x 3 ,ya)} 

- "12 Kbyb'( x ^y^ x ^yi> T cdJd'( x 2,y2; x 2 ,y 2 )r ef J f ,{x 3 , y 3 ; x s , y 3 Wa'b',c'd',e' f ,{x 1 ,y 1 ;x 2 ,y 2 ; x 3 , y 3 ) . 

(60) 

I 



There are two types of objects that we need to compute 
in order to use this relation: cumulants of the micro- 
scopic densities that are generated by the differentiation 
of In Z m with respect to v a and v a b , and second deriva- 
tives of T m with respect to p a and h ab . In the end, we 
want to evaluate these objects in the liquid phase where 
h is equal to zero. But we know that the free energy 
can be written as the HNC free energy plus 2PI contri- 
butions, that are 0(/i 3 ), i.e. contain more than three h 
links. Thus when taking one or two derivatives of the 
2PI diagrams with respect to h, they still will contain at 
least one h link, and will all cancel out when evaluated 
at zero h. This proves that Eq. (60) when evaluated in 
the liquid phase can be computed by replacing the T m 
functionals in the r.h.s. by r^ NC . 

Thus the needed derivatives of T m can be computed 
starting from Eq. (34) by dropping the 2PI diagrams. 
The only difficult term is the sum of ring diagrams. Its 
first derivative is given in Eq. (38). This equation can 
be rewritten with the propagators and vertex functions 
G and rW, and can then be easily differentiated with 
respect to h by noting that, since they are inverse from 
each other, we must have: 



r (2) 



1 

G 



5IH 



5G 



1 

CP 



_ r (2)p(2) 



(61) 



Explicitly, we obtain 

sr<$(x uyi ) _ 

SG c d(x 2 ,y 2 ) 
Sc a b(xi,yi) 



rZ\ x uyiKd'( x 2,y2), 

S ab 6(x 1 ,y 1 ) - c afc (xi,yi) 



^(2), 



Sh cd (x 2 ,y 2 ) \p a (xi) 
1 



I — -, — rS c dS{x 2 ,y 2 ) - c cd {x 2 ,y 2 ) 

Pc(X 2 ) 

(62) 

One can also obtain this result by a diagrammatical anal- 
ysis, since c is written as a sum of chains of h bonds. We 
thus obtain: 



T a)cd( X i^ X 2,y2) \ Hq = -\p{8acS{x U X 2 ) + 5 ad 5{x U y 2 )) , 

(63) 

^ a °b%( x ^y^ x 2,y2) \ Uq 

= ^P 5 ac 5 M V\^{xi,x 2 )T\^{ y i, y2 ) (64) 

= ^P 5 ac 5 bd r^ c (a;i,X2;yi,2/2) • 

We see that these derivatives are delta functions for what 
concerns the replica indexes. If we drop the space in- 
dexes, we can thus already perform trace over replica 
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indexes in Eq. (60), to get: to obtain: 

p(0,3) = ]_ W (3,0) 
ab,cd,ef ^3 a,c,e 

^rft ® wgj, - {two perms.} 

+ -/hNC ® ^HNC ® + ^ tW ° P CmS -> 

-8r (0 ' 2) «r (0 ' 2) «r (0 ' 2) «w (0 ' 3) 

(65) 

where ® means space convolution with respect to two 
spatial indexes. 

Finally, we only have left the task to compute the 
derivatives of the logarithm of the partition function (the 
terms 

ty(o,3) j w {1 ' 2 \...). These terms, equivalently to 
the one written in Eq. (18) for the liquid, are averages 
of microscopic one or two-point densities, that are easily 
computed within the 1RSB framework that we use here, 
as we explain in the following. 

B. 1RSB structure of the theory 

We denote by (•) the equilibrium average for the repli- 
cated system. Once again, in the limit v a b(r) — > 0, all 
replicas fall in the same state but are otherwise uncorre- 
cted inside the state. Finally, we want to evaluate all our 
averages in the liquid phase. This leads to the following 
rule to compute the average (•): one should 

• factorize the averages (•) when they involve differ- 
ent replicas, and 



• remove the replica indexes. 

• replace (•) = (») liq 

For instance, for any spatial argument, and for a ^ 6, we 
have that following the prescription above 

(PaPb) = (Pa) (Pb) = (P) liq (p) Hq = P 2 • (66) 

Similarly, assuming that different letters denote different 
values of the indexes: 

(PaPaPb) = (PaPa) (Pb) = (PP) Hq (p) Uq = P{Ghq + p 2 ) ■ 

(67) 

PabPac) = (PaPb PaPc) = {PaPa) (M (Pc) 

= {p)u q (p)u q (PP)u q ( 68 ) 
= p 2 {G Hq + P 2 ) . 

We will thus obtain quantities that do not depend on 
replica indices anymore, allowing to sum over these in- 
dices, and finally take the m — > 1 limit. The free-energy 
will have an overall factor m(m — 1), and thus we will 
consider the free energy divided by m(m — 1). All the 
terms that are 0((m — l) 2 ) or higher-order in to — 1 will 
thus disappear in the limit m — > 1. This will eliminate a 
large number of terms in the complicated expression (60). 
Indeed, recalling that we will calculate a free-energy cor- 
rection of the form: 



5T m [p,h] = / F i °fJ def (x 1 ,y 1 ;x2,y2;x 3 ,y 3 )h ab (x 1 ,y 1 )h cd (x 2 ,y2)hef{x 3 ,y 3 ) 

I r ( 69 ) 

EEE 



Cy^d Cy^f 



^ab%,ef( x tiy^ x zi 2/2; x 3 , y 3 )h{x 1 ,y 1 )h(x 2 , y2)h{x 3 ,y 3 ) 



xi,yi,--- ,x 3 ,y 3 



r 



we see that the trace over replica indexes will let appear 
to (for example the number of ways to choose the index 
a) and to — 1 factors (the number of ways to choose b =/= a, 
once a has been chosen). If the terms in r^ ' 3 * 1 did not 
depend on replica indexes at all, we would obtain a factor 
(in — l) 3 from the three summations over replica indexes 
6, d, and /, and the contribution would cancel out in 
the to — y 1 limit. In the different terms of Eq. (65), we 
see that the first two lines contain quantities that do not 
depend on two or three indexes. These terms will give 
contributions that are at least 0((m— l) 2 ), and thus will 
cancel out in the m — > 1 limit. Another remark allows 



to discard also the terms coming from the third line of 
Eq. (65): each time a replica indexes structure for W^ 1 ' 2 ' 
is chosen, it is easy to see that this structure always have 
at least a multiplicity to(to— 1), thus making these terms 
C((to - l) 2 ) also. 



Finally we are left with only the last line of Eq. (65) to 
compute, which amounts to compute W(°' 3 \ under the 
prescription described in the beginning of this subsection. 
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C. Final calculation 

Before turning to the explicit evaluation of W^-°' 3 \ it 
is useful to remark that, within the 1RSB structure that 
we have, we can parametrize the replica dependance of 
our quantities in a simple way. 

Take a matrix that depends on two pairs of replica in- 
dexes M ab<cd ,, with a ^ b and c ^ d (we do not explicit 
the space indexes). Examination of the different possi- 
bilities for a, b, c and d shows that we have only three 
genuinely different possibilities: 



M, 



ab,cd 



^ab,ab 



(70) 



This is a consequence both of the 1RSB ansatz and of the 
symmetry of the functions with respect to permutations 
of indexes (and their associated space indexes) . This can 
be summarized in: 



M, 



b.cd 



--Mi 



, , Sac + Sad + he + 5bd . , T 
M 2 V 1VL 3 



(71) 



where Mi , M 2 and M 3 are related to the above terms by: 

l/i = 2 [M oM & - 2M oMc + M o6 , cd ] , 
M 2 = 4[M ab , ac - M ab , cd ] , (72) 
M 3 = M ab , cd . 

I 



The quantity that we are interested in is a matrix that 
depends on three pairs of indexes. In this case there are 
8 topologically different possibilities [32] : 



W. 



(0,3) 



ab,cd,ef 



W\ = W a b,ab,ab , 
W 2 = W a b,bc,oa , 
W 3 = W ab ,ab,ac , 
Wa = W a b,ab,cd , 
W 5 = W ab ,ac,ad , 
W 6 = W ab ^ ac . bd , 
W 7 = W a b,ac,de , 
W 8 = W a b,cd,ef ■ 



(73) 



A relation like Eq. (71) is again possible, but cumber- 
some, and we do not write it explicitly because we do 
not need it. By using the prescription for calculating av- 
erages of one and two-point densities described above, we 
can easily compute the Wi. We obtain: 



Wi(xi,--- ,2/3) 



W 2 (xt,-- 
W 3 (x u .- 



■ »J/3) 



W$ (x 1 ,x 2 ,x 3 ) W$ (2/1,2/2,2/3) 

+pW l f a ) (x 1 ,x 2 ,x 3 ) (p 2 + Gu q (y 1 ,y 2 ) + G Hq {yi,y 3 ) + G Hq (y 2 ,y 3 )) 

+pW^ q (yi,y2,y 3 ) (p 2 +Gu q (x 1 ,x 2 ) + G Uq (x 1 ,x 3 ) + G Uq (x 2 ,x 3 )) 

+P 2 (Gu q (x 1 ,x 2 )Gi iq (yi,y 3 ) + Gii q {x 1 ,x 2 )Gu q (y 2 ,y 3 ) + Gu q (xi, x 3 )Gu q {yi, y 2 )) 

+p 2 (Gu q (x 1 ,x 3 )Gu q (y 2 ,y 3 ) + Gu q (x 2 ,x 3 )G Hq (yx,y 2 ) + Gi iq (x 2 ,x 3 )Gu q (y 1 ,y 3 )) 

G Uq (xi,x 2 )G Uq (?/i , x 3 )G Uq (1/2 , 2/3 ) 

+P 2 (Guqixx^^Gugiyi^xs) + Gi iq {xi,x 2 )G Hq {y 2l y 3 ) + Gu q (yi,x 3 )G Hq {y 2 ,y 3 )) 
P 3 W[*}(x 1 ,x 2 ,x 3 ) + pGu q (yi,y 2 ) (w ( - 3) (xi,x 2 ,x 3 ) + pG Uq {x u x 3 ) + pGu q (x 2 ,x 3 ) 



, (74) 



• ,Vz) = ^P 2 G ilg (xi,a; 2 )G ;ig (2/i,a;3) 

• ,Z/s) = ^P 3 W l f q ) {x 1 ,x 2 ,x 3 ) , 



(75) 



and we find We = W7 = W% = 0. The factors 1/8 come 
from the fact that a derivative of In Z m with respect to 
— (3v gives p^ /2 and not p( 2 \ Of course, the choice 
of spatial indexes is arbitrary, and one can make any 
permutation, as long as it respects the symmetry of the 
functions. We can now perform the trace over replica 



indexes, which will give, apart from the multiplicity of 
each of the Wi, additional factors due to the m 4 1 
limit: a term which has a factor m — 2 will have an extra 
— 1 factor, a term that has a factor (m — 2)(m — 3) will 
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have an extra factor 2, etc.. We obtain: 

^ M^hr) E E E W K,e/ 

v ; a^6 c^d e^/ 

= AWi - 8W 2 - 24VF 3 + 16VF 4 + 16W 5 



(76) 



W l ^(x 1 ,x 2 ,x 3 )W^(y u y 2 ,y 3 ) 
-2Giig(x 1 ,x 2 )Gii q (y 1 ,x 3 )Giiq{y2,y3) 



Finally we can perform the convolution with the deriva- 
tives of the HNC free energy in Eq. (65), and by making 
repeated use of the second- and third-order OZ equations 
we obtain: 



ST m [p,h] 

= -£ 
6 



V(x 1: --- ,y3)h(xi,yi)h(x2, y 2 )h(x 3 , y 3 ) , 

xi,--- ,3/3 



V( Xl , 



>Vs) = ^u ) q {xi,x 2 ,x 3 )T^ ) q {y 1 ,y 2 ,y 3 ) 

- r S Oi » ^a)r|2 (yi , a?s)rS2 (2/2 , 2/3 ) 



(77) 



This third order correction includes the RHNC term, 
which reads (here V is the system's volume): 



p 2 V 
12 



hirf 



p 3 V 



h(q) 3 



q (1 + P%)) 3 



(78) 



Substracting it from the result in Eq. (77), we obtain 
finally: 



m(m — 1) 



^HNC 



[pM 



m(m — 1) 



EL 
12 

6 



c (3) (xi,x 2 ,x 3 )c (3) (2/i,2/2,2/3)Mxi, 2/i)Mx2, 2/2)^(^3, 2/3) 



(79) 



c (3) (xi , x 2 , x 3 )h(x 1 , y)h(x 2 ,y)h(x 3 ,y) . 



Xl,X2,X3,y 



r 



This is the desired result: the next order term in the 
order-parameter expansion beyond the RHNC approxi- 
mation. From this approximation of the free-energy, we 
can make use of the variational principle Eq. (33) to ob- 
tain a closed equation on h: 

c(r) = h(r) - ln[l + h(r)} 

c Uq ( r ' f i' r 3)c (3) (0, r 2, r4,)h(ri, r 2 )h(r 3 , r 4 ) 



£ ,(3) 



c^>{r,r 1 ,r 2 )h(r 1 )h(r 2 ) 



2 ' 



c (3) (0,ri,r 2 )/i(r - n)h(r - r 2 ) 



(80) 



which provides the first correction to Eq. (43). Using 
translational invariance as well as the invariance under 
permutation of the three variables of we get: 



~(3)/ 



(r-i,r 2 ,r 3 ) = c (3) (ri - r 2 ;ri - r 3 ) 



-(3), 



r; s) 



with r = ri — r 2 s = r\ — r 3 



„(3) 



(r;s) 



(s,r) 



„(3) 



(— r; s — 7') 



Defining the double fourier transform of as: 



c^(k,q) = I e 

> r,s 



'(r,«) , 



we obtain two invariance principles: 

c (3) (fc, g) = c (3 ) (4, fc) = c^ 3) {-k-q,q) 



(81) 

(82) 
(83) 



Performing a Fourier transformation on our equation we 
get 

c(fc) = T (h - ln[l + h] \ (fc) 
+ tt [ c (*)(-k,q)cW(k,-q)h(q)h(k-q) 



c^(q,k-q)h(q)h(k-q) 



c^\-q,-k + q)h{q)h(k-q) , 



(84) 



which using the invariances in Eq. (83) is simplified to: 
c{k) =F(h- ln[l + h]) (k) 



l + p 2 cS 3 \q,k-q)\ -l)h(q)h(k-q) 



(85) 

Before turning to a numerical resolution of this equation, 
we can again make the naive expansion of the HNC term 
in powers of h, and keep only the lowest order term, to 
obtain: 

= SB [ p ^(k,-qfS( q )S(k-q)f(q)f(k-q) . 



l-f(k) 2p J q 

(86) 

This recovers exactly the three-body term in the MCT 
kernel Eq. (48). 
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V. THREE-BODY CORRELATIONS 
AND NUMERICAL SOLVING 

We thus have obtained a closed equation of the order- 
parameter, that necessitates as an input the two- and 
three-body direct correlation functions of the liquid. We 
already quoted that we decided here to work within 
the PY approximation for the two-point functions. It 
is known [24] that the PY approximation [39], which 
amounts to treat the fluid in a mean-field approximation, 
but under the exact constraint that the pair correlation 
function should vanish for distances smaller than 1 [46], 
is particularly efficient for hard spheres. Furthermore, we 
dispose of an analytic expression for the two-body direct 
correlation function in that approximation [37, 38] . How- 
ever, the three-body direct correlation function still needs 
to be approximated. Computing the third-order direct 
correlation function is in itself a hard problem of liquid 
theory. The best approximation available was shown by 
numerical works [34] to be the HNC3 approximation de- 
veloped by Attard [33]. However, this approximation is 
very computationally demanding, and because our pur- 
pose here is merely to demonstrate the importance of 
higher-order terms in the expansion in powers of the or- 
der parameter, we do not aim at quantitative efficiency, 
and wish to find a simpler approximation scheme. 



A. Denton &: Ashcroft approximation 

A good compromise between simplicity and efficiency 
for evaluating the third-order direct correlation func- 
tion [35] is the Denton- Ashcroft approximation [36] . This 
approximation gives an analytic form that necessitates as 



an input only the second-order direct correlation func- 
tion, for which we can use the PY result. Within their 
approximation, c' 3 ) is given by: 



.0) 
"DA 



(k,q) 



1 



c^(k)d p c^(q) + c^(q)d p c^\k) 



C (2)(0) 

_ V 2) (Q) ( 2)(fc)c ( 2)( )- 

(c( 2 )(0)) 2 



(87) 



Within this approximation, angular dependance is ne- 
glected. We can recover it by symmetrizing the expres- 
sion: 



.(3) 
"DAS 



(k,q) 



.(3) 



DA (k,q) + c ( §' A (k, \k + q\) 



+cf A {q,\k + q\) 



As stated before, we use as an input the Percus-Yevick 
direct correlation function, that reads (in units of the 
hard-sphere diameter): 



cf Y {r) 



1 o 

6(fbr — —(par r < 1 , 



(89) 



% r > 1 . 

where (p is the packing fraction defined by <p = 7rp/6, 
a = (1 + 2<p)7(l - iff and b = -(1 + <p/2) 2 /(l - <p) 4 . 
The corresponding Fourier transforms, and derivatives 
with respect to density are simply computed analytically. 
Using this approximation for c*- 3 -* , and writing the inte- 
grals in bipolar coordinates by using the isotropy of the 
liquid, we obtain the following set of equations: 



m 

F(r) 
H(k) 



F(k) + H(k) , 
h(r) - In f 1 + h(r) 



pS{k) 
8n 2 k 



du 



k+u 



\k—u\ 



dv u vh(u)h(v) 



(90) 



P 2 c { D AS {u,v]k) 



- 1 



where we have defined 



CDAs( u > v 'i k ) = ^ 



c%\(k, u) + c% A (k, v) + c% A (u, v) 



(91) 

Equations (87) - (91) now completely specify our ap- 
proximation. We present in the following a preliminary 
numerical resolution in order to demonstrate the impor- 
tance of the correction H . 



B. Numerical resolution methodology 

The usual calculation for HNC alone is more stable if 
we write the iteration procedure in terms of c and x = 
h — c, which would give in our case: 



c(r) 



x(r)+H(r) 



1 - X(r) , 



(92) 



or the same equation with H = in the case of RHNC. 
We first solve the RHNC equation Eq. (90) with H = 0, 
by using PY approximation for the two-point functions. 
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without 3 body term 
with 3 body term 




without 3 -body term 
with 3 -body term 




FIG. 4: Non-ergodicity factor as a function of the wave vector 
very close to the critical point tp — ipa, for the PY + HNC 
calculation, and with inclusion of the three body term. 



FIG. 5: Non-ergodicity factor as a function of the wave vector 
at fixed packing fraction 0.6, for the PY + RHNC calculation, 
and with inclusion of the three body term. 



Once a stable solution of the RHNC equations has been 
found we introduce the three-body correction H, and 
solve by the very same Picard iterative scheme that is 
used for solving HNC and RHNC. Explicitly, the resolu- 
tion will give: 

• Start from the value of kc(k), obtained with the 
previous iteration, or from the old HNC value if 
first iteration 

• Use Eq. (40) to deduce the value of kh{k) and k\(k) 

• Use kh(k) to evaluate kH(k) 

• Inverse Fourier transform k\(k) and kH{k) to ob- 
tain r\(r) and rH(r) 

• Use it to evaluate the new rc(r) with Eq. (92) 

• Fourier transform rc(r) to obtain kc(k) 

• Mix it with the old kc(k) to avoid rapid changes 

• Repeat these steps until fcc(fc) has converged 

We used a grid of 2 10 equally spaced points on a box of 
size 11, and a mixing parameter 0.01. We note that eval- 
uation of the correcting term H in Eq. (90) has a compu- 
tational cost of order of the square number of points on 
the grid, significantly slowing down the resolution of the 
equation, since to avoid instabilities, c is made to evolve 
very slowly by the mixing procedure. 

C. Results and discussion 

We solved our improved equation (90) for several pack- 
ing fractions starting from <p = 0.6, and decreasing it un- 
til the non-trivial solution for h disappeared. We found 
that the inclusion of the three-body terms, even when us- 
ing a crude approximation such as the symmetrized Den- 
ton & Ashcroft approximation, leads to a strong shift of 



the dynamical transition point from <pd = 0.591897 to 
tfd = 0.555860, which is closer to the predicted MCT 
transition but farther from the numerically estimated 
transition. It is hard to guess what would be the evo- 
lution of ipd when additional corrections are added. 

The resulting non-ergodicity factors at the transition 
are depicted in Fig. 4. The inclusion of the three body 
term has significantly enhanced the first peak of f(k) 
upon inclusion of our correction. These preliminary 
results must be treated with caution, because it was 
found [42] that the numerical solving of the RHNC equa- 
tion is sensitive to the discretization used, and very large 
grids with a very large number of points must be used in 
order to obtain stable results. This situation is expected 
to be the same with the presence of the 3 body term, 
but the 0(N 2 ) scaling of the numerical resolution in our 
case prevented us from performing a stability analysis. 
The qualitative picture is nevertheless not expected to 
be modified by these considerations. 

In Fig. 5 we show the results inside the glass phase, 
at a fixed packing fraction (p = 0.6, with and without 
inclusion of the three-body term. In this case the effect 
of including the 3-body term is much bigger and goes in 
the right direction of increasing f(k) at small k (although 
still not enough at very small k below the first peak). 

VI. HIGHER ORDERS 

Our calculation provides the exact free-energy and cor- 
relation function at order h 2 . This feature allows us to 
put into correspondence the MCT kernel, which is also 
C(/ 2 ), with the replica result. We have seen that the 
MCT kernel in Eq. (48) contains exactly the three-body 
contribution we obtain with replica theory, which is al- 
ready a quite surprising result, given the differences that 
exist between the two approaches. However, it shows that 
the two-body terms contained in the MCT kernel are thus 
impossible to obtain within a static framework. The pc- 
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culiar wave-vector dcpcndance of these terms arise from 
the calculation of forces, inherent in dynamical theories, 
but absent in static ones. One could however wonder 
whether these two-body terms could arise in a dynami- 
cal calculation because of the factorization of a dynami- 
cal four-point vertex function, that would be forced to be 
expressed as a C(/ 2 ) function. Indeed, the main approx- 
imation involved in Mode-Coupling theories is the factor- 
ization of a four-point function, and since MCT breaks 
down at high dimensions, it is possible that in the pro- 
cess, "glassy" correlations are factorized along with "liq- 
uid" ones, forcefully introducing these new C(/ 2 ) terms. 
This is nevertheless highly speculative, and no satisfy- 
ing dynamical theory exists yet, that would be able to 
investigate these considerations. 

An interesting feature of our calculation is that, once 
the theory is set up, we can already uncover the next 
terms with a diagrammatical visualization of the expan- 
sion. To go further, we can now wonder what is the 
term 0(h 4 ) in the free-energy (that would correspond to 
a 0(f 3 ) term for the correlation function). The require- 
ment that we must have 2PI diagrams is quite strong, 
and it is easy to convince one-self that the only possi- 
ble diagram that we can construct is shown in Figure 6. 
The four-point functions that are connected by h lines 
are made of 1PI diagrams, and we will thus obtain a new 
contribution to the free-energy that contains the 4-point 
vertex functions of the liquid: 

const x/ rW(ri,...,r r )r[5(rs I ... > r 8 ) 

Jr 1 ,...,r 8 (93) 

x h(n,r 2 ) ■ ■ - Hrr^s) . 

Interestingly, we don't even have to work out the precise 
diagrammatics behind this procedure, since RHNC gives 
a contribution in every diagram, and at all orders in h, 
so that we can use RHNC to fix the prefactor of each 
diagram we compute. 




FIG. 6: Diagrams that contribute to the free-energy at order 
h 4 . A wiggly line joining two replica indices a and b is a 
h a b, with a ^ b function, a black dot attached to a zone with 
replica index a is an integration point weighted by a density 
factor p a - 

Computing the n-th order direct correlation function is 



still a difficult problem, but with approximations such as 
the Denton- Ashcroft approximation, that are issued from 
density functional theory, we can deduce from the ap- 
proximation for c*- 3 -* the corresponding approximation for 
c(") by successive differentiation with respect to density. 
The last difficulty is that an 0(f n ) kernel will require the 
numerical evaluation of an n-dimensional integral, that 
have a computational cost of order 0(N n ). 

In order to visualize the increasing difficulty of going 
to the next orders, we show in Figure 7 the possible dia- 
grams at order h 5 . We see that new, intricate terms arise, 
and that in general, the 0(h n ) contribution will contain 
all vertex functions of the liquid of orders ranging from 
3 to n. 




FIG. 7: Diagrams that contribute to the free-energy at order 
h 5 . A wiggly line joining two replica indices a and b is a 
h a b, with a ^ b function, a black dot attached to a zone with 
replica index a is an integration point weighted by a density 
factor p a - 



VII. CONCLUSION AND DISCUSSION 

Our analysis shows that the 2PI corrections to the 
RHNC free energy are quantitatively relevant at the dy- 
namical transition, and must be properly taken into ac- 
count in order to obtain an accurate static description 
of the transition. It seems then that only by doing this 
properly we will be able to make a clear connection be- 
tween dynamical and static theories of the dynamical 
transition of glasses. In this respect our results are the 
following: 
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• At the level of the RHNC approximation, the equa- 
tion for f(k) can be developed in powers of f(k). 
The result is Eq. (47), which corresponds to the 
Mode-Coupling equation with a kernel equal to 1. 
The latter is also the result of a zero-th order field 
theory calculation as reported in [9] . 

• The first correction to the RHNC approximation 
provides an additional contribution to the kernel 
which happens to correspond exactly to the three- 
body term of the Mode-Coupling kernel. 

• The next corrections will give terms proportional 
to / 3 in the right hand side of Eq. (47) , hence no 
additional contributions to the Mode-Coupling ker- 
nel can be generated by these terms. We are forced 
to conclude that there is no way of generating the 
terms of the Mode-Coupling kernel proportional to 
c(k) by means of a static computation. 

• We have shown the way to compute higher-order 
terms, although we expect the numerical resolution 
of the corresponding approximations to be hard. 

It would be therefore very important to perform a similar 
calculation (namely, a systematic expansion in powers 
of /(&)) also on the dynamical side. This would allow 
for a systematic comparison of the results. One would 
then obtain a proper theory for the crgodicity breaking 
that occurs at pj, free of the ambiguities of MCT, and 
systematically improvable. Work has been done in this 
direction in the last years [9, 11, 47, 48], but the situation 
is still unsatisfactory. 

1000 f 



without 3-body term 
with 3-body term 




0.1 



0.01 



0.1 
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On the static side, the first task to be performed is to 
reconcile the small cage expansion and the 2PI approach 
that we use in this work in order to have a unified static 
theory. The way to do this is indicated by the shape of 
the pair correlation function between two different repli- 
cas at the transition, depicted in Fig. 8. We see that g 
develops a strong peak at r = 0, that dominates the rest 
of its features. This peak simply reflects the fact that, in 
the glass phase, two replicas tend to be very close from 
each other. In the small cage expansion, this idea is re- 
flected in the introduction of a cage size parameter A, 
which serves as a expansion parameter. In our case a 
direct introduction of the cage size A would be more dif- 
ficult, because it would amount to parametrize the pair 
correlation function by the cage size, which is difficult to 
do without specifying a given shape of cages. In the high 
dimensional limit, the cage can safely be approximated 
by a Gaussian as far as the free energy is concerned [49] , 
but this is not expected to hold in our low dimensional 
case [50]. 

It is interesting to note that we have observed, as in 
Fig. 8 that the inclusion of the three-body term system- 
atically increases the spatial separation between the first, 
glassy, peak of g(r) and the subsequent, liquid peaks. A 
clear separation between these contributions allows for 
an unambiguous definition of the molecules introduced 
in order to perform small-cage expansions [13, 15], and 
is clearly a sign that reconciling RHNC with small-cage 
expansions may not be out of reach. 

The large contribution coming from the peak at r = 
in g shows that the diagrams that contribute the most 
to the free-energy are the most connected ones [15]. It 
should be possible to put these diagrams in correspon- 
dence with the diagrams re-summed in the small cage 
expansions of [13] and [15] in order to make progress. 



FIG. 8: Inter-replica pair correlation function g(r) at packing 
fraction ip — 0.6, with and without three-body correction. 
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